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Highly accurate results for the homogeneous electron gas (HEG) have only been achieved to date 
within a diffusion Monte Carlo (DMC) framework. Here, we introduce a newly developed stochastic 
technique, Full Configuration Interaction Quantum Monte Carlo (FCIQMC), which samples the 
exact wavefunction expanded in plane wave Slater determinants. Despite the introduction of a basis 
set incompleteness error, we obtain a finite-basis energy which is significantly, and variationally lower 
than any previously published work for the 54-electron HEG at r s — 0.5 a.u., in a Hilbert space of 
10 108 Slater determinants. At this value of r s , as well as of 1.0 a.u., we remove the remaining basis 
set incompleteness error by extrapolation, yielding results comparable or better than state-of-the-art 
DMC backflow energies. In doing so, we demonstrate that it is possible to yield highly accurate 
results with the FCIQMC method in sizable periodic systems. 

PACS numbers: 71.10.Ca,31.15.V-, 71.10.-w, 71.15.-m 



The homogeneous electron gas (HEG), described by 
a simple model Hamiltonian, encapsulates many of the 
difficulties with modern electronic structure theory. To 
date the only truly successful ah initio methods to yield 
accurate ground state energies at a range of densities 
have been quantum Monte Carlo techniques, in partic- 
ular diffusion Monte Carlo (DMC) [IHf|. The most fa- 
mous of these was the results of Ceperley and Alder from 
which the LDA functionals of Density Functional Theory 
were parameterised [l|, 0] ■ Diffusion Monte Carlo would 
be an exact technique but for the fixed-node approxima- 
tion, which requires the nodes in the wavefunction due 
to fermionic exchange to be specified in advance by some 
trial wavefunction. In general, the fixed-node approxima- 
tion lacks a method of being systematically improved to 
find the exact result. Attempts to go beyond the fixed- 
node approximation have been met with some success, 
however complete elimination of this error has not been 
achieved [H 0, 0, @]. In particular, the release node (RN) 
method is practical only in systems for which the Bosonic 
ground state is close in energy to the Fermionic one. In 
the HEG this is only true at low density. At high den- 
sities, the RN-DMC is unstable, and fixed-node DMC 
with backflow corrections is the most viable option. This 
leaves open the question of the magnitude of the remain- 
ing fixed-node error. 

Full configuration interaction (FCI) aims to find the 
wavefunction expressed as a linear combination of Slater 
determinants, formed from rearrangements of N elec- 
trons in an underlying one-electron basis of M spatial 
orbitalsQ, Q. This is equivalent to an exact diagonaliza- 
tion of this space. Since such a basis set of Slater deter- 



minants scales as 



M 
.N/2) 



benchmarks from FCI are ex- 



tremely difficult to produce. There has been surprisingly 
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little work undertaken with polynomially-scaling high- 
accuracy quantum chemical techniques, even though it 
has recently been shown that finite systems ranging from 
as few as 54 electrons can begin to capture the physics of 
the 3D HEG accurately 1(1 U|- In part this might be due 
to the required size of the one-electron basis and that, on 
approaching the thermodynamic limit for metals, many 
approximate methods find divergent energies [12|. In con- 
trast, truncated configuration interaction will tend to- 
wards zero correlation energy. 

We present the application of a new method, FCI quan- 
tum Monte Carlo (la-Hit, which stochastically samples 
the exact wavefunction providing the accuracy of ex- 
act diagonalization at a greatly reduced computational 
cost, to the high-density 54-electron HEG at r s =0.5 and 
1.0 a.u. This is the regime in which backflow corrections 
to FN-DMC are the largest [1 % . 

The Model- We seek to find the ground-state wave- 
function and total energy of the iV-electron HEG 
simulation-cell Hamiltonian: 



H 



-Nv M 



(1) 



where the two-electron operator v a p is the Ewald inter- 
action, 



0, q=0 1 ' 



i>m is the Madelung term, which represents contribu- 
tions to the one-particle energy from interactions between 
a point charge and its own images and a neutralising 
background 10> 12 , 18 1 , and Q is the real-space unit cell 
volume. 

We use an expansion of the wavefunction in a Slater 
determinant basis, 
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(3) 
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where each determinant is a normalised, antisymmetrizcd 
product of plane waves, 



A = A [V ! i(Xi)^(x 2 )...'(/'fc(xAr)] 



(4) 



(5) 



The i index, which uniq uely labels each determinant, is 
its normal-ordered string 19]. The wave vectors k are cho- 
sen to correspond to the reciprocal lattice vectors of a 
real-space cubic cell of length L, 



(6) 



where n,m and I are integers. The Hartree-Fock determi- 
nant is the determinant occupying N plane waves with 
the lowest kinetic energy. The full basis set for our calcu- 
lation is constructed of all Slater determinants that can 
be made from M plane waves (2M spin orbitals) forming 
a closed-shell of orbitals in k-space with a kinetic energy 
lower than an energy cutoff ^k^. Plane waves are conve- 
nient because taking a single cutoff parameter to infinity 
makes the one-electron basis set complete. Moreover, 
plane waves are natural orbitals for the electron gas, im- 
plying that a FCI expansion is rapidly convergent in this 
basis [2d]. 

The determinant expansion given in Eq. [3] can be 
inserted into the imaginary-time Schrodinger equation, 
yielding a set of coupled equations for the determinant 
coefficients 



dd 
dr 



(Ha - S)Ci 



(7) 



Setting dC-Jdr — and solving for S by exact diagonal- 
ization yields the total energy for the problem in a given 
basis. 

In a novel and recently developed quantum Monte 
Carlo, termed Full Configuration Interaction QMC 
(FCIQMC)[I3], Eq.[7]is regarded as a set of master equa- 
tions governing the dynamics of the evolution of the de- 
terminant coefficients in imaginary time, with elements 
of H being non-unitary transition rates. These dynamics 
are simulated by introducing a population of N w 'walk- 
ers' distributed over the determinants, which arc signed 
to represent the sign of the coefficients within the simu- 
lation, C\ oc (N\ (t) } . The walker population is then al- 
lowed to evolve through discretized imaginary time-steps 
by spawning, death/cloning and annihilation events ac- 
cording to Eq. [7] until a steady-state is reached. The 
exact rules for this can be found in [l3| . 

The parameter S, termed the shift, is a population con- 
trol parameter which can be updated self-consistently at 
equilibrium to oscillate around the total energy. How- 
ever, throughout this work, the projected energy is used 
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(a) A typical i-FCIQMC run. At r ~ 3.8 a.u., the shift S was 
allowed to vary to keep the walker number at an average of 20 
million. From this point, an average was taken of the total energy. 
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(b) i-FCIQMC calculations with n a< j<j=3 are run at increasing N w 
values, with the aim that the limit N w — > oo is found by the 
simulation at maximum walker number. 

FIG. 1: Plots showing calculation of the 
i-FCIQMC energy for N = 54, 2M = 682, r s = 0.5 a.u. 
The result is reached in the limit of long time (iteration 
number) and large walker number N w . 
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M S ir.l stochastic correlation energy estimator, 

j 



Sfciqmc = (E(t)) = £ (A |ff I A)) 



(No(t))' 



(8) 



where Do is taken as the Hartree-Fock determinant and 
the sum j need only be taken over the O [7V 2 M] doubly 
excited determinants of Do • 

Typically the system is initially grown by setting S to 
some positive value and allowing evolution from a single 
determinant to allow an unbiased evolution of the popula- 
tion. Only populations above a critical system-dependent 
size are able to converge to the FCI distribution, and this 
size scales linearly with the size of the Hilbert space [13 1. 

In order to alleviate this problem, an adaptation of this 
method has been developed, called initiator-FCIQMC 
(z-FCIQMC)[li- 16 1 . The determinant space is instanta- 
neously divided into those determinants exceeding a pop- 
ulation of n a dd walkers, termed initiator determinants, 
and those that do not. When considering a determinant 
whose current population is zero, the sum in the second 
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FIG. 2: Convergence with walker number up to 
approximately 100 million walkers (taken to be the 
infinite walker limit) is shown for a variety of basis sets 
for N = 54, r s = 1.0 a.u. Each line is labelled with the 

spinorbital number, 2M and was calculated with 
n a dd=3. As the basis set size grows, so the size of space 
and the number of walkers required to sample the space 
accurately grows. 



term of Eq. the term describing net flux of walkers 
onto that determinant, is taken to be only over initiator 
determinants. This effectively introduces a survival of 
the fittest criterion for survival of newly spawned walk- 
ers. If a walker has been spawned from a determinant 
with an instantaneous population exceeding a parame- 
ter n a( jd, the child is allowed to survive. However, if 
the parent walker is a determinant with a population 
smaller or equal to n a dd then the child only survives if it 
has been spawned to a currently occupied determinant. 
This i-FCIQMC has been shown to dramatically acceler- 
ate the convergence of FCIQMC with respect to walker 
number. Note that in the large walker number limit, the 
i-FCIQMC tends to the FCIQMC algorithm, which it- 
self converges rigourously to the FCI energy. Figure Q] 
illustrates an i-FCIQMC energy calculation in this way. 
Previous work has shown that the rapid convergence to 
this limit can be examined by finding correlation energies 
at increasing walker numbers (Fig. Ilbp fl6|. 

As the basis set size grows, so the number of walkers 
required to recover the total energy to a given level of 
accuracy increases (Fig. Further results are taken to 
be converged with respect to this initiator error. 

Basis set extrapolation.- Although i-FCIQMC is able 
to produce exact results in a finite basis set, these are only 
upper-bound estimates to the true ground-state energy. 
This error in the energy, termed the basis set incomplete- 
ness error, is absent in DMC results which do not have a 





FCI 2M = 2838 


FCI M = oo 


1 a.u. 


/ a.u. per electron 


/ a.u. per electron 


0.5 


3.22086(2) 


3.2202(2) 


1.0 


0.53073(4) 


0.5300(3) 



TABLE I: i-FCIQMC total energies for 2M spin 
orbitals. The error estimate for the finite-basis 
corresponds to stochastic error. The M = oo result is 
based on extrapolations shown in Fig. [3J from which 
the error estimate derives. 



substantial dependence on a basis set [21 j. 

Extrapolation of the correlation energy to the com- 
plete basis set limit is performed regularly in molecu- 
lar systems for which scaling laws have been investi- 
gated extensively [22j]. In plane wave systems, a 1/M 
extrapolation is used for the basis set incompleteness er- 
ror in methods employing the random phase approxima- 
tion and second-order M0ller-Plesset theory 23J, |24| . We 
note that analytic expressions can be derived with these 
methods for the HEG that also show a 1/M relation- 
ship. Figure [3] illustrates that by using this fit at high 
basis set sizes, complete basis set exact diagonalization 
energies can be obtained that compare well with most 
recent high-accuracy DMC results [6J. The unquantificd 
initiator error is sub-mEh, and is the largest source of 
error in these results but these are nonetheless thought 
to be upper-bound estimates of the exact energy, due to 
the observed variationality of the initiator error in large 
basis sets (see Fig. [2]). 

Results and Conclusions.- Results of 
i-FCIQMC calculations performed on the 54-clcctron 
gas, for basis sets containing between 162 and 2838 
spinorbtials, are shown in Fig. |3] and Table HI In these 
calculations, Hilbert spaces ranging from 10 39 to 10 108 
Slater determinants (Fig. [3]) were sampled to produce 
high-accuracy energies for the HEG using approximately 
300,000 core hours in total. 

For r s =1.0 a.u., we obtain a variational finite-basis 
result that lies just below the fixed-node DMC result, 
with the extrapolated energy agreeing well with back- 
flow DMC energies. At this r s , the remaining fixed node 
error can be estimated by variance extrapolation of the 
backflow DMC results, and is thought to be sub-mEhQ. 
The results presented here, containing a similar order of 
magnitude of error, seem to substantiate this claim. To 
resolve the remaining fixed-node error in the backflow re- 
sults, we would need to further reduce the leading source 
of error in the i-FCIQMC results, which is due to basis 
set incompleteness. 

For r s =0.5 a.u., we obtain a variational finite-basis 
result that lies below the backflow DMC result. How- 
ever, our extrapolated energy falls significantly below 
the lowest DMC result found to date, suggesting resid- 
ual fixed-node error in the backflow DMC energies is of 
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FIG. 3: z-FCIQMC total energies for a basis of 2M spin orbitals. Each basis set corresponds to a kinetic energy 
cutoff, with 2M = 2838 corresponding to 208 Ryd at r s ~0.5 a.u. and 52.1 Ryd at r s =1.0 a.u. Each calculation used 
40 million walkers for r s = 0.5 a.u. and 100 million walkers for r s = 1.0 a.u.. The blue dashed line is an 
extrapolation to M — > oo based on the expected form 1/M using the data set with the largest number of walkers, 
shown with error bars in the inset. The DMC results, taken from Rfos et a/.0], do not suffer from basis set error and 
are shown as two horizontal lines representing the mean plus and minus one standard deviation. Almost identical 
backflow results can be found for r s = 1.0 a.u. in a study by Kwon et al.\^. 



the same order as the backflow corrections themselves. 
It has been commented that the wavefunctions for the 
3D HEG change much more at higher densities than at 
lower densities under backflow transformations Our 
results suggest that in spite of this greater change, the 
backflow transformation still comparatively struggles at 
r s = 0.5 a.u. 

The significance of these results extends beyond the 
sheer size of the many-electron basis which is being effec- 
tively sampled without error. The fact the results can be 
taken as exact within the designated basis set allows them 
to be used to benchmark other, more approximate meth- 
ods in this system, as well as providing upper bounds to 
which the DMC method can optimise its nodal surface. 
Our method can also easily be extended to examine other 
properties of the HEG, in particular momentum distri- 
butions and Fermi liquidparameters, which is the focus 
of many current studies [2h6|. Ill|. Beyond the HEG, the 
results also bode well for the treatment of more realistic 
solid-state systems, including metals, at unprecedented 
accuracy. 

Acknowledgements.- The authors gratefully acknowl- 
edge Neil Drummond for many useful conversations and 
to Mike Towler for providing us the CASINO code for 
comparison benchmark values at the early stages of this 
studv(25l|. This work was supported by a grant from the 
Distributed European Infrastructure for Supcrcomputing 
Applications under their Extreme Computing Initiative. 



[1] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 

(1980) . 

[2] M. Holzmann, D. M. Ceperley, C. Pierleoni, and K. Esler, 

Phys. Rev. E 68, 046707 (2003). 
[3] G. Ortiz and P. Ballone, Phys. Rev. B 50, 1391 (1994). 
[4] Y. Kwon, D. M. Ceperley, and R. M. Martin, Phys. Rev. 

B 58, 6800 (1998). 
[5] I. G. Gurtubay, R. Gaudoin, and J. M. Pitarke, J. Phys.: 

Condens. Matter 22, 065501 (2010). 
[6] P. Lopez Rfos, A. Ma, N. D. Drummond, M. D. Towler, 

and R. J. Needs, Physical Review E 74, 066701 (2006). 
[7] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 

(1981) . 

[8] J. Olsen, B. Roos, P. j0rgensen, and H. Jensen, J. Chem. 

Phys. 89, 2185 (1988). 
[9] P. J. Knowles and N. C. Handy, Chem. Phys. Lett. Ill, 

315 (1984). 

[10] N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. 

Foulkes, Physical Review B 78, 125106 (2008). 
[11] M. Holzmann, B. Bernu, C. Pierleoni, J. McMinis, D. M. 

Ceperley, V. Olevano, and L. Delle Site, Phys. Rev. Lett. 

107, 110402 (2011). 
[12] L. Onsager, L. Mittag, and M. J. Stephen, Annalen der 

Physik 473, 71 (1966). 
[13] G. H. Booth, A. J. W. Thorn, and A. Alavi, J. Chem. 

Phys. 131 (2009). 
[14] D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. 

132, 041103 (2010). 
[15] D. Cleland, G. Booth, and A. Alavi, J. Chem. Phys. 134, 



5 



024112 (2011). 

[16] G. H. Booth, D. Cleland, A. J. W. Thorn, and A. Alavi, 

J. Chem. Phys. 135, 084104 (2011). 
[17] P. Ewald, Ann. Phys. 64 (1921). 

[18] L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. 

Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev. B 

53, 1814 (1996). 
[19] W. Kutzelnigg and D. Mukherjee, J. Chem. Phys. 107, 

432 (1997). 

[20] E. Davidson, Rev. Mod. Phys 44, 451 (1972). 

[21] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra- 



jagopal, Rev. Mod. Phys. 73, 33 (2001). 
[22] W. Kutzelnigg and J. D. Morgan, J. Chem. Phys 96, 
4484 (1992). 

[23] J. Harl and G. Kresse, Physical Review B 77, 045136 
(2008). 

[24] A. Griineis, M. Marsman, and G. Kresse, J. Chem. Phys. 

133, 074107 (2010). 
[25] R. J. Needs, M. D. Towler, N. D. Drummond, and P. L. 

Rfos, J. Phys.: Condens. Matter 22, 023201 (2010). 



